Solution for (1+1) dimensional surface solitons in thermal nonlinear media 
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Analytical solutions for (1+1)D surface fundamental solitons in thermal nonlinear media are 
obtained. The stationary position and the critical power of surface solitons are obtained using this 
analytical solutions. The analytical solutions are verified by numerical simulations. The solutions 
for surface breathers and their breathing period, and solutions for surface dipole and tripole solitons 
are also given. 
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I. INTRODUCTION 

In nonlinear optics, an optical beam forms a soliton 
when the self-focusing caused by the nonlinearity of the 
medium balances the diffraction. Recently, people pay 
more attention to nonlocal nonlinear media 'l'-'4'| . There 
are some types of solitons in nonlocal media, such as 
multipole-mode solitons [^-01^^^^®'^'^^ Hermite soli- 
tons and vortex solitons 0, 

Surface waves have been studied generally in physics, 
chemistry and biology. They have been used to study the 
surface properties of media and the interaction between 
media and optical beams. In the presence of nonlinearity, 
some kinds of surface solitons have been found theoret- 
ically and experimentally, such as in local Kerr media 
jl2| , waveguide arrays Il3l4l5l , photorefractive media 
Uq . and metamaterials 19] 
Nowadays, surface solitons in nonlocal nonlinear me- 
dia have been found theoretically and experimentally 
[20| . An interesting property of surface solitons is that a 
beam launched away from the stationary position moves 
to the interface and oscillates at its vicinity, even if the 
launch position is far away from the stationary position. 
Now, varies types of surface solitons have been found and 
studied, such as incoherent surface solitons [2l|, surface 
dipoles [12, and surface vortices [12] . Surface dipoles 
are found to be stable in the entire existence domain. 
Surface vortices are found to exhibit strongly asymmet- 
ric intensity and phase distributions. Bound states of 
surface vortex solitons belong to a class of surface soli- 
tons having no counterparts in bulk media [25] . 

However, there are no analytical solutions for the fam- 
ily of surface solitons due to the boundary conditions. 
We notice that the optical field of surface solitons almost 
resides in the nonlinear medium if the refractive index 
difference between two media is large enough. It is rea- 
sonable to assume that the optical intensity equals zero 
at the interface with large refractive index difference. In 
this paper, we obtain analytical solutions for the family 
of surface solitons and surface breathers under this as- 
sumption. Using the analytical solutions, we also obtain 



the stationary position and the critical power of surface 
solitons, and the breathing period of surface breathers. 
These results are proved by numerical simulations. 



II. SURFACE FUNDAMENTAL SOLITONS 

First, we consider a (l-l-l)D sample sketched in Fig. 
[IJa) following the experiment in Ref. [1^. The sam- 
ple width d is 1 mm. The interface {x = d) between 
a nonlocal nonlinear medium (lead glass) with a refrac- 
tive index ni = 1.8 and a linear medium (air) with a 
refractive index n2 = 1.0 is thermally insulating, thus 
the temperature distribution T satisfies the boundary 
condition dT / dx\x=d = 0. The other side of the sam- 
ple is thermally conductive at the temperature Tq, Tq is 
the temperature in the absence of light, thus the bound- 
ary condition is T\x=o = Tq. The nonlinear refractive 
index change (An) of the sample satisfies the relation 
An = l3AT = /3(T - To), where /3 is the thermal coeffi- 
cient of the refractive index. An satisfies the conditions 
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T satisfies the heat diffusion equation [3| 



(1) 



(2) 



where n is the thermal conductivity coefficient, / is the 
optical intensity, a is the absorption coefficient. 

We consider the propagation of a TE polarized laser 
beam along z axis in the vicinity of the interface. In the 
nonlocal nonlinear medium (0 < x < d), the slowly vary- 
ing light field amplitude U{x,z) (i.e. / — \U\^) satisfies 
the nonlocal nonlinear Schrodinger equation. 
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where ki — riiko, kg is wave number in vacuum and 
|An| <ti ni. In the linear medium (x > d), U{x,z) satis- 
fies 
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where ^2 = n2ko. Here we ignore the absorption of the 
medium, which is necessary for producing temperature 
changes and thermal nonlinearity. If the absorption is 
considered in propagation, the rigorous unchanged soli- 
ton does not exist in the nonlocal nonlinear medium due 
to the loss of the beam energy and the decreasing of the 
nonlinearity [23, H^- Fortunately the propagation dis- 
tances in most experiments are short enough so that the 
influence of the absorption can be ignorable d, [2^ [2l| . 
In order to show the stability of solitons in this paper, 
the propagation distances in numerical simulations are 
chosen to be 100-500 diffraction lengths, which are much 
longer than that in real experiments. 

In Fig. [TJa), the solid line is the optical intensity pro- 
file for the beam launched at x = xq and the dashed line 
is the nonlinear refractive index which is not a bell shape. 
Obviously, the optical beam will approach to the inter- 
face where the nonlinear refractive index is larger, which 
presents the nonlinear attraction of the interface. The 
optical beam will be reflected when it arrives at the in- 
terface because the linear refractive index of the sample 
is larger than that of the air. When the nonlinear at- 
traction balances the interface reflection, and simultane- 
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FIG. 1. (Color online) (a) Sketch of a (1+1)D surface soliton. 
(b) Propagation of a surface soliton, and (c) its intensity and 
nonlinear refractive index profiles, (d) Propagation of a first- 
order Hermite-Gaussian soliton in a bulk medium, and (e) its 
intensity and nonlinear refractive index profiles. The beam 
width is 30 /j,m and Z_r — kiWQ is the Rayleigh distance. The 
widths of the media are 1 mm for the surface soliton, and 2 
mm for the bulk soliton, respectively. The solid line is the 
optical intensity profile and the dashed line is the nonlinear 
refractive index profile. 



ously the nonlinear self-focusing balances the diffraction, 
the beam propagates stably along z axis, which forms 
a surface soliton as shown in Figs. [Ijb) and[IJc). Be- 
cause the refractive index difference between two media 
is large enough, i.e. ni — n2 = 0.8, there is very little 
optical intensity residing in the air. In the simulation, 
the ratio of the optical intensity at the interface {x = d) 
to the maximum optical intensity is approximately 10~^, 
so the optical intensity at the interface can be negligible. 
Therefore we introduce the boundary conditions 



0, n 



0, 



(5) 



where I\x=o = is valid because the optical beam is nar- 
row and far away from the left boundary. Using above 
conditions and Eq. ([T]) for the temperature distribution, 
one can completely solve Eqs. ^ and ^ for the prop- 
agation of optical beams in the nonlinear medium, i.e. 
< a; < d. Here the physical influence of the linear 
medium, governed by Eq. Q with continuous conditions 
at the interface, is simplified to the boundary condition 
I\x=d — for the circumstance when the refractive index 
difference between two media is large enough. 

To obtain analytical solutions for Eqs. ([2]) and (O, 
we consider a (1+1)D first-order Hermite-Gaussian (HG) 
soliton in a same bulk medium with the sample width 2d, 
as shown in Figs. [Hd) and[lje). The optical intensity 
satisfies the boundary conditions 
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0, 
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because the optical beam is narrow and far away from 
the two boundaries. It also satisfies the condition 



I\:r=d = 0, 



(7) 



because the first-order HG function is antisymmetric and 
launched at the sample center {x = d). The temperature 
T satisfies the conditions T\x=o = Tq and T\x=2d = Tq, 
because the two boundaries of the sample are thermally 
conductive at the temperature Tq. Therefore the refrac- 
tive index change satisfies the boundary conditions 
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The refractive index change is symmetrical for a symmet- 
rical optical intensity, so it satisfies 
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The evolution of the first-order HG soliton in the bulk 
medium (0 < a; < 2d) is governed by Eqs. ([2|) and ([3]) 
under the boundary conditions Eqs. (j6|) and ([8]). Com- 
paring Eqs. ([7]) and (|9]) for the antisymmetric HG bulk 
soliton with Eqs. ([T|) and ([5]) for the surface soliton, one 
can find that the surface soliton is identical with the half 
part (0 < a; < d) of the first-order HG bulk soliton. 
Therefore, the solution of the surface soliton can be ob- 
tained from the solution of the first-order HG soliton in 
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the bulk medium, namely, 
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where Hi is the first-order Hermite polynomial, wq is 
the beam width, b is the propagation constant, Uq is a 
constant, and < x < d. 
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FIG. 2. (Color online) The intensity profile of surface solitons 
for (a) Wq — 15 and (b) wo = 60 fim. The square symbols 
represent the numerical solutions and the solid lines represent 
the analytical solutions, (c) and (d), propagations of surface 
solitons using the analytical solution as incident profile corre- 
sponding to (a) and (b), respectively. The sample width is 1 
mm. The dashed lines indicate the interface. 



Regarding the analytical solution as a trial solution, we 
will find the numerical soliton profile using the iterative 
method based on Eqs. In our simulations, the pa- 

rameters are a = 0.01 cm^^, Tq = 25 °C, k = 6.37 x 10"^ 
WK~^cm~\ /3 = 1.4 X 10~^ K~\ A = 514 nm, which are 
the same as those in Ref . ^2^] . The numerical profiles are 
calculated without the assumption I\x=d = 0, whereas 
all results of intensities at a; = d are small enough to be 
approximated to zero. Figs. [2][a) and[2jb) show the com- 
parison of the optical intensity of the analytical solutions 
(solid lines) with the numerical solutions (square sym- 
bols) for Wo — 15 ^m [Fig. [UJa)] and wo = 60 ^m [Fig. 
Iljb)]. One can see that the analytical solutions are in 
good agreement with the numerical solutions. Figs. HJc) 
and IHd) show the propagations of surface solitons us- 
ing the analytical solutions Eq. (ITOl) as incident profiles. 
The propagations are stable for a relative long distance, 
which verifies our analytical solutions again. 

The critical power of surface solitons can be obtained 
from the analysis of bulk solitons. In the bulk medium, 
the nonlinear refractive index in Eq. ^ can be expanded 
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where Pq is the input power of bulk solitons. The pa- 
rameter (5^ can be found to be proportional to the beam 
intensity at a; = d, or be proportional to Po/[\/ttw{z)]^ 
i.e. 



a/3 



0rw(z) 

Then, Eq. (|3]) reduces to a linear equation 
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By substituting the solution of bulk solitons [same as 
Eq. (Uni) but < a; < 2d ] into Eq. (US]), and from the 
coefficient of each order term of x, we obtain the relation 
between the soliton power and the beam width in the 
bulk medium, 



Pc 



(14) 



The critical power of surface solitons is half of that of 
bulk solitons, i.e. Pg = Pc/2. The amplitude constant 
t^o = VPc/ {y^woY^^ , and the propagation constant 



(15) 
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The stationary position Xs, namely the maximum op- 
tical intensity position of surface solitons, can be also 
obtained. 



Xs ^ d — Wq. 



(16) 



Now the analytical solution for fundamental surface soli- 
tons is obtained completely. For example, in Fig. ^c), 
the beam width wq is 15 /im, the critical power P^ is 
888.4 W/cm [in the simulation P^ = 892 W/cm] and the 
stationary position is 15 fim away from the interface. In 
Fig. Hfd), the beam width wo is 60 /im, the critical power 
Ps is 13.9 W/cm [in the simulation Ps = 14.2 W/cm] and 
the stationary position is 60 /im away from the interface. 

Figure[3ja) shows the relation between the beam width 
and the critical power of surface solitons. The solid line 
is the theoretical value calculated by Eq. (|14p and the 
square symbols are the numerical value. It can be seen 
that the theoretical value is in good agreement with the 
numerical value. The relation between the beam width 
and the propagation constant b is shown in Fig. [3l^b). 
One can see that there exists a great difference between 
the analytical solution and the numerical results. It 
means that Eq. ^T^i is incorrect because we neglect 
the higher order terms of the expansion of the nonlin- 
ear refractive index in Eq. (1111) . The higher order terms 
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on Eqs. ©-(HI). The input condition is Ut = U + MJ = 
U\\ + ei?(a;)], where A[/ is the perturbation amphtude, 
Ut is the total amphtude, U is given by Eq. (ITUl) . e 
is a perturbation constant, which is 5% in our simula- 
tion. R{x) is a random function whose value is between 
and 1. Figure Hl^a) shows the initial intensity distribu- 
tion of the surface soliton with a great distortion. After 
propagating about 1 Rayleigh distance as shown in Fig. 
m^b), the distortion disappears and the intensity profile 
becomes smooth. Figure Ul^c) shows the intensity distri- 
bution after propagating 100 Rayleigh distances. This 
intensity distribution is still smooth and without distor- 
tion. Compared with Fig. HJ^b), there is little change of 
the optical intensity profile in Fig. IDJc). So, the distorted 
profile can reform itself into a proper soliton after a short 
distance (1 Rayleigh distance) and propagate along the 
sample surface stably. 



FIG. 3. (Color online) (a) Relation between the critical power 
Ps and the beam width wo. (b) Relation between the beam 
width wq and the propagation constant b. (c) Relation be- 
tween the soliton critical power Ps and the propagation con- 
stant b. The solid line is the theoretical results and the square 
symbols are the numerical results. 



have little influence on the beam intensity and the crit- 
ical power, but a great deal of influence on the phase 
change and the propagation constant. Because of the 
same reason, the dependence of the critical power Pg on 
the propagation constant h does not agree with the nu- 
merical solution, as shown in Fig. El^c). 
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FIG. 4. (a) The initial intensity distribution with 5% noise. 
The intensity distributions of (b) after propagating 1 Rayleigh 
distance and (c) after propagating 100 Rayleigh distances. 
The beam width is 30 /xm. 

In order to confirm the stability of our solutions, we 
simulate the beam propagation with perturbation based 






FIG. 5. The sectional profile of the optical intensity of the 
surface soliton propagation for different sample widths, (a) 
The sample width d = 0.5 mm (~ IGwo) , (b) d = 0.09 mm 
{3wo), and (c) d = 0.06 mm {2wo). The beam width is 30 
fj,m. 

Finally, we discuss the influence of the sample width 
on the surface solitons. Fig. [TJb) shows the propaga- 
tion of the surface soliton for the beam width wq — 30 
/im and the sample width d = 1 mm (~ 33wo), which 
is far greater than the beam width. If the sample width 
is reduced to d ^ 0.5 mm idwo) as shown in Fig. 
OJa), the propagation of the surface soliton is not influ- 
enced by the left boundary. And if the sample width 
is reduced to d — 0.09 mm which is three times of the 
soliton beam width, the influence of the left boundary is 
still small. However when the sample width is reduced 
to d = 0.06 mm {2wo) as shown in Fig. [SJc), the prop- 
agation of the surface soliton is unstable and there is a 
great influence of the left boundary on the optical beam, 
and the optical beam can not form a stable soliton. We 
take the three times of beam width as a critical value. 
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When d/wQ = 3, the optical intensity at left boundary is 
about I\x=o/ Imax ~ 3 X 10"'^, and the influence of the 
left boundary is small enough to be negligible. When 
the sample width is smaller than three times of the beam 
width, I\x=o is not negligible, e.g. I\x=Q/Imax ~ 0.1 for 
d/wo — 2, and surface soliton can not exist due to the 
influence of the left boundary. This property is identical 
with that of solitons in bulk media 261. 



III. SURFACE BREATHERS 

An optical beam can form a breather when the input 
beam power is stronger or weaker than the soliton critical 
power in the bulk medium. Fig. [Hb) (Pq = 2Pc) and Fig. 
ini^d) (Pq = -Pc/2) show the first-order HG breathers with 
Wo = 30 /im in bulk media. Since the antisymmetry of 
the first-order HG beams in bulk media can retain during 
propagation, it is reasonable to expect that an optical 
beam launched at the sample surface also owns breathing 
effect. When the surface beam power is stronger than 
the soliton critical power, i.e. Pq = 2Ps as shown in 
Fig. mja), the optical beam approaches to the interface at 
first and converges because of the limit of the boundary. 
As the beam arrives at a position near the interface, it 
is reflected by the interface and diverges again. When 
the beam is far away from the interface, it is attracted 
again to the interface. This property is identical with 
the first-order HG breather in the bulk medium shown 
in Fig. m^b) at < a; < d. If the surface beam power is 
weaker than the soliton power, i.e. Pq = Ps /2 as shown 
in Fig. EJc), the breathing effect is opposite to the case 
of the stronger power. The optical beam is reflected by 
the interface and diverges at flrst, then attracted by the 
interface and returns back. 

The analytical solution of surface breathers can be also 
obtained from the result of bulk breathers. In the bulk 
medium, we search for a breather solution for Eq. Q in 
the HG function form 



U{x) = t/o-ffi 



[x-d) 
w{z) 



[x - df 
2w{zf 



exp — — I + ic{z){x — d)^ 

where c(z) is the wavefront curvature. Substituting the 
breather solution into Eq. ([T^. and from the coefficient 
of each order term of x, we obtain two equations ^27,] 



dw{z) 2c{z)w{z) 
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The combination of Eq. (jl8b[) with the derivative form 
of Eq. yields 



1 



[w{z)/wo] 

d(z/Zfl)2 2[w{z)/wo]^ 



^=0, (19) 
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FIG. 6. (Color online) The propagations of surface breathers 
and bulk breathers, (a)and (c), surface breathers with Pq = 
2Ps and Po = 0.5Ps, respectively, (b)and (d), HG breathers 
with Po = 2Pc and Pq = 0.5Pc, respectively. The initial beam 
width is wo — SOfim for all. (e) The relation between the 
surface breathing period and Po/Ps for the analytical value 
(solid line) and the numerical value (square symbols). 



Using the method introduced in Rcf. 



we obtain 



w{z) = Wo 



VA- 



5A(A- 1) 




2AZr 



(20) 

where A = (2Pc/Po)^' Then the parameters c{z) and b 
can be gotten and the solution of HG breathers in bulk 
media is obtained. The analytical solution of surface 
breathers is just the half part of Eq. ([T7)) . i.e. < a; < d. 
It should be noted that the critical power of surface 
breathers is a half of that of bulk solitons, i.e. Pg = Pc/2 



From Eq. 
tained as 



20]), the breathing period Lb can be ob 
2V6ttZii /2P,^ 



L 



Pn 



(21) 



A surface breather has a same period as its bulk coun- 
terpart. In Fig. inije), the analytical relation (solid line) 
approximates well to the numerical results (square sym- 
bols). 
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IV. SURFACE MULTIPOLE SOLITONS 
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FIG. 7. (Color online) (a) The profile of a dipole surface 
soliton. The square symbols are the numerical solution and 
the solid line is the analytical solution, (b) Propagation of 
the dipole soliton using the analytical solution as the incident 
profile. (c)Propagation of the dipole soliton using the numer- 
ical solution as the incident profile. The sample width is 1 
mm and the beam width is 250 fim. 

In this section, we discuss the surface multipole soli- 
tons using the same method mentioned above. A surface 
dipole soliton can be expressed as a half part {0 < x < d) 
of a third-order HG soliton in the center of a bulk medium 
(0 < X < 2d). The profile of a dipole soliton (solid line) is 
shown in Fig. [Tj^a), and the numerical iterative solution 
(square symbols) is also shown for comparing. It can be 
seen that there exists a little difference between the nu- 
merical solution and the analytical solution. It is known 
that surface dipole solitons are stable [23i], analogously 
quadropole solitons in bulk media are stable too [l^ . 
The difference between the numerical solution and the 
approximate analytical solution can be considered as a 
perturbation, which leads to some spatial oscillations in 
propagation when the analytical solution is used as the 
initial profile, as shown in Fig. [Tljb). However, the prop- 
agation using the numerical solution is stable as shown 
in Fig. Ulc). 

A surface tripole soliton can be expressed as a half part 
(0 < a; < d) of a fifth-order HG soliton in the center of 
a bulk medium (0 < x < 2d). The profile of a tripole 
soliton (solid line) is shown in Fig. [DJa), and the numer- 
ical iterative solution (square symbols) is also shown for 
comparing. There also exists a little difference between 
the numerical solution and the analytical solution. It 
has been proved in Ref. that surface tripole solitons 
are unstable. As shown in FigEl^c), the oscillation oc- 
curs after propagating over relatively long distance, i.e. 
z/Zf( w 100, when we use the numerical solution as the 



incident profile. When we use the approximate analyti- 
cal solution, the perturbation resulted from the difference 
between the numerical solution and the approximate an- 
alytical solution grows up explosively. As a result, the 
strong oscillation takes place during the propagation as 
shown in Fig. [8][b). Anyway, our analytical solution co- 
incides approximately with the numerical solution and it 
will be helpful for analyzing the multipole surface soliton 
propagation. 
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FIG. 8. (Color online) (a) The profile of a tripole surface 
soliton. The square symbols are the numerical solution and 
the solid line is the analytical solution, (b) Propagation of 
the tripole soliton using the analytical solution as the inci- 
dent profile. (c)Propagation of the tripole soliton using the 
numerical solution as the incident profile. The sample width 
is 1 mm and the beam width is 120 /im. 



V. CONCLUSION 

In conclusion, the analytical solution for surface funda- 
mental solitons has been found by comparing with bulk 
solitons, also for surface breathers, surface dipole and 
tripole solitons. The key assumption in this paper is 
that the optical fields all locate in the nonlinear medium 
when two media have large difference in the linear re- 
fractive index, and the intensity at the interface can be 
negligible as J = 0. This method is similar to the image 
beam method [28j in dealing with the boundary force 
exerted on spatial solitons, which is an analogue of the 
well-known image charge method. By introducing image 
beams, the boundary conditions are automatically satis- 
fied and analytical solutions can be found conveniently. 
The numerical simulations show our assumption is rea- 
sonable and correct. This method can be used also in 
dealing with other surface solitons if the linear refractive 
index difference between two media is large enough. 
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